FAIRE DES CARTOGRAMS DANS R

Lorem Ipsum is simply dummy text of the printing and typesetting industry. Lorem Ipsum has been the industry’s standard dummy text ever since the 1500s, when an unknown printer took a galley of type and scrambled it to make a type specimen book. It has survived not only five centuries, but also the leap into electronic typesetting, remaining essentially unchanged. It was popularised in the 1960s with the release of Letraset sheets containing Lorem Ipsum passages, and more recently with desktop publishing software like Aldus PageMaker including versions of Lorem Ipsum.

Packages

installer les packages suivants ….. + version min

install.packages(knitr)
install.packages(sf)
install.packages(mapsf)
install.packages(packcircles)
install.packages(cartogram)
install.packages(recmap)
install.packages("https://cran.r-project.org/src/contrib/Archive/cartogramR/cartogramR_1.0-1.tar.gz", repos = NULL, type = "source")

Import et mise en forme des données

données INSEE (pop + CSP en 2018)

Le package sf bla bla bla…..

library(sf)
communes <- st_read("data/isere.geojson", quiet = TRUE ) %>% st_transform(2154)
data <- read.csv("data/popisrere.csv",  dec = ",")
communes = merge(x = communes[,c("id","name","geometry")], 
                 y = data[,c("id", "pop2018","agri", "art", "cadr", "interm", "emp","ouvr","retr")],
                 by = "id")
isere = st_union(communes)
knitr::kable(communes[c(0:10),], row.names = F, digits = 1)
id name pop2018 agri art cadr interm emp ouvr retr geometry
38001 Les Abrets en Dauphiné 6336 30.1 210.9 256.0 712.9 833.3 948.8 1397.2 MULTIPOLYGON (((900716.2 64…
38002 Les Adrets 1025 0.0 53.0 168.7 163.9 82.0 67.5 159.1 MULTIPOLYGON (((931354.5 64…
38003 Agnin 1127 0.0 28.9 72.2 168.6 139.7 134.9 192.6 MULTIPOLYGON (((844459.1 64…
38004 L’Albenc 1279 25.0 30.0 65.0 215.0 130.0 175.0 210.0 MULTIPOLYGON (((893395.8 64…
38005 Allemond 954 5.0 75.0 30.0 160.0 150.0 120.0 170.0 MULTIPOLYGON (((934142.3 64…
38006 Allevard 4062 0.0 176.7 323.2 535.2 469.6 469.6 984.2 MULTIPOLYGON (((948125.2 64…
38008 Ambel 28 10.0 0.0 0.0 0.0 0.0 0.0 10.0 MULTIPOLYGON (((930843.9 64…
38009 Anjou 1009 10.0 54.9 60.0 124.9 100.0 115.0 278.5 MULTIPOLYGON (((847739.9 64…
38010 Annoisin-Chatelans 686 14.8 44.4 54.2 113.3 64.1 78.9 113.3 MULTIPOLYGON (((875494.4 65…
38011 Anthon 1073 5.0 55.0 75.0 185.0 125.0 130.0 220.0 MULTIPOLYGON (((866538.7 65…

mapsf

Le package mapsf permet de faire des cartes thématiques dans R. C’est le package qui succède au package cartography.

Création d’un template cartographique

library(mapsf)

col = "#c291bc"
credits = paste0("Bronner Anne-Christine & Nicolas Lambert, 2021\n",
                  "Source: IGN & INSEE, 2021")
theme = mf_theme(x = "default", bg = "#f0f0f0", tab = FALSE, 
                   pos = "center", line = 2, inner = FALSE, 
                   fg = col, mar = c(0,0, 2, 0),cex = 1.9)

template = function(title, file, basemap = TRUE, scale = TRUE){

  mf_export(
    communes,
    export = "png",
    width = 1000,
    filename = file,
    res = 96,
    theme = theme, 
    expandBB = c(-.02,0,-.02,0)
)
  
if (basemap == TRUE){
  mf_shadow(x = communes, col = "grey50", cex = 1, add = TRUE)
  mf_map(communes, col ="#CCCCCC", border = "white", lwd = 0.5, add = TRUE)
}

mf_title(title)

if (scale == TRUE){
  mf_scale(size = 20, pos = "bottomright", lwd = 1.2, cex = 1, col = "#383838", unit = "km")
}

  mf_credits(
    txt = credits,
    pos = "bottomleft",
    col = "#1a2640",
    cex = 0.8,
    font = 3,
    bg = "#ffffff30"
  )
}
template("Template cartographique", "maps/template.png")
mf_map(communes, col = col, border = "white", lwd = 0.5, add = TRUE)
dev.off()

Symboles proportionnels

template("Symboles proportionnels (mapsf)", "maps/prop.png")
mf_map(communes, var = "pop2018", col = col, border = "#6b4266", type = "prop",
       inches = 0.8, leg_title_cex = 1.2, leg_val_cex   = 0.8,
       leg_title = "Nombre d'habitants, 2018")
dev.off()

template("Symboles proportionnels (mapsf)", "maps/prop2.png")
mf_map(communes, var = "pop2018", col = col, border = "#6b4266", type = "prop",
       inches = 0.8, leg_title_cex = 1.2, leg_val_cex   = 0.8, symbol = "square",
       leg_title = "Nombre d'habitants, 2018")
dev.off()

Dot density maps

dotdensitymap <- function(x, var, onedot = 1, radius = 1){
x <- x[,c("id",var,"geometry")]
x[,"v"] <- round(x[,var] %>% st_drop_geometry() /onedot,0)
dots <- st_sample(x, x$v, type = "random", exact = TRUE)
circles <- st_buffer(dots, dist = radius)
return (circles)
}
onedot = 500
dots = dotdensitymap(x = communes, var = "pop2018", onedot = onedot, radius = 300)
template(paste0("Un point = ",onedot," habitants"), "maps/dotdensity.png")
mf_map(dots, col = col, border = "#520a2c", lwd = 0.5, add = TRUE)
dev.off()

packcircles

Le package packcirles propose 3 algorithmes simples pour déplacer des diques sur un plan 2D de telle sorte qu’ils ne se supperposent pas. Nous pouvons l’utiliser pour créer des cartogrammes de Dorling (Dorling 1996).

Création d’un ficher de données simplifié avec les coordonnées des centroides des communes.

dots = communes
st_geometry(dots) <- st_centroid(sf::st_geometry(dots),of_largest_polygon = TRUE)
dots <- data.frame(dots$id, dots["pop2018"], st_coordinates(dots))
dots = dots[,c("dots.id","X","Y","pop2018")]
colnames(dots) <- c("id","x","y","v")
dots <- dots[!is.na(dots$v),]
knitr::kable(dots[c(0:5),], row.names = F, digits = 1)
id x y v
38001 902026.9 6495943 6336
38002 934555.0 6466630 1025
38003 845708.9 6472468 1127
38004 892892.6 6460697 1279
38005 937029.6 6456880 954

La fonction circleRepelLayout() prend un ensemble de cercles dans un cadre de données et utilise la répulsion itérative pour essayer de trouver un arrangement sans chevauchement où tous les centres des cercles se trouvent à l’intérieur d’un rectangle de délimitation. Si aucun arrangement de ce type ne peut être trouvé dans le nombre maximum d’itérations spécifié, la dernière tentative est renvoyée.

library("packcircles")

k = 500 # pour ajuster la taille des cercles
itermax = 10 # nombre d'iterations

dat.init <- dots[,c("x","y","v")]
dat.init$v <- sqrt(dat.init$v * k)
simulation <- circleRepelLayout(x = dat.init, xysizecols = 1:3,
                                wrap = FALSE, sizetype = "radius",
                                maxiter = itermax, weights =1)$layout
knitr::kable(simulation[c(0:5),], row.names = F, digits = 1)
x y radius
902120.7 6496140 1779.9
934555.0 6466630 715.9
845708.9 6472468 750.7
892892.6 6460697 799.7
937029.6 6456880 690.7
circles <- st_buffer(sf::st_as_sf(simulation, coords =c('x', 'y'),
                      crs = sf::st_crs(communes)), dist = simulation$radius)
circles$v = dots$v

template("Dorling (packcircles)", "maps/dorling1.png", scale = FALSE)
mf_map(isere, col = "#CCCCCC", border = "white", lwd = 0.5, add = TRUE)
mf_map(circles,col = col, border = "#6b4266", lwd = 0.5, add = TRUE)
dev.off()
## png 
##   2

cartogram

**le package cartogram est développé par Sebastian Jeworutzki. Il propose trois méthodes : Dorling, Olson et Dougenik

library(cartogram)

Dorling

dorling = cartogram_dorling(communes, "pop2018", k = 1.8)
template("Dorling (cartogram)", "maps/dorling2.png", scale = FALSE)
mf_map(isere, col = "#CCCCCC", border = "white", lwd = 0.5, add = TRUE)
mf_map(dorling, col = col, border = "#6b4266", lwd = 0.5, add = TRUE)
dev.off()

Olson

Bla bla bla cartogrammes discontinus…. (Olson 1976)

x : un objet sf weight : nom de la variable k : facteur pour augmenter la taille des polygons inplace : Si VRAI, chaque polygone est modifié à sa place initiale, si FAUX, les multi-polygones sont centrés sur leur centroïde initial

olsen <- cartogram_ncont(communes, "pop2018", k = 1, inplace = TRUE)
template("Olsen (cartogram)", "maps/olsen.png", basemap = FALSE, scale = FALSE)
mf_map(isere, col = "#CCCCCC30", border = "white", lwd = 0.5, add = TRUE)
mf_map(olsen, col = col, border = "#6b4266", lwd = 0.5, add = TRUE)
dev.off()

Douguenik

Bla bla bla…. (Dougenik, Chrisman, and Niemeyer 1985)

x : un objet sf weight : om de la variable itermax : nombre d’itérations maximum si maxSizeError n’est pas atteint maxSizeError : la déformation s’arrete si l’erreur moyenne est inférieur à cette valeur prepare : mettre “adjust” permer d’acceler le temps de calcul. threshold : seuil pour la préparation des données

Dougenik <- cartogram_cont(communes, "pop2018", prepare = "none", itermax = 10, maxSizeError = 1.15)

Calcul des érreurs

sumarea = sum(as.numeric(st_area(Dougenik)))
sumpop = sum(Dougenik$pop2018)
Dougenik$error = (as.numeric(st_area(Dougenik)) /  sumarea) / (Dougenik$pop2018 / sumpop) * 100
summary(Dougenik$error)
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
##    0.2513   96.9756  109.4084  143.6154  130.7821 2734.5953
bks = c(min(Dougenik$error),70,80,90,100,110,120,max(Dougenik$error))
cols = c("#d53e4f", "#f46d43","#fdae61","#fee08b","#e6f598","#abdda4", "#66c2a5")

Affichage de la carte

template("Dougenik (cartogram)", "maps/dougenik.png", basemap = FALSE, scale = FALSE)
mf_map(x = Dougenik, type = "choro",var = "error", pal = cols, breaks = bks, border = "#6b4266", lwd = 0.5, add = TRUE)
dev.off()

RecMap

Recmap (Heilmann et al. 2004) (Panse 2016) est un package développé par Christian Panse. Il perlet de concertir les géométries des unités spatiales en rectangles dont la surace est définie en fonction d’une donnée statitique sauntitative absolue. L’algorithe de RecMap est disponible ici. Développé en C++11, le package dépend des packages GA (>= 3.1), Rcpp (>= 1.0), sp(>= 1.3)

Attention, il ne fonctionne que sur un nombre limité d’unités territoriales. Et même dans ce cas, il ne fonctionne pas très bien, notamment avec les petites unités spatiales.

library(recmap)

Ici un exemple sur les pays d’Europe.

europe <- st_read("data/europe.geojson")

Préparation des données

coords = data.frame(st_coordinates(st_centroid(st_geometry(europe))))
bb <- lapply(st_geometry(europe), function(x){st_bbox(x)})
dx <- unlist(lapply(bb, function(x){x[3]-x[1]})) / 2
dy <- unlist(lapply(bb, function(x){x[4]-x[2]})) / 2

df <- data.frame(x = coords$X, 
                 y = coords$Y, 
                 dx = dx, 
                 dy = dy, 
                 z = europe$pop2008,
                 name = europe$id)
knitr::kable(df[c(0:10),], row.names = F, digits = 1)
x y dx dy z name
4631849 2726828 281676.5 144612.8 8318592 AT
3947159 3073426 130414.4 108636.1 10666866 BE
5560592 2309244 241662.7 187871.0 7640238 BG
4190230 2634505 175597.7 112221.5 7593494 CH
6424605 1654151 89911.7 81559.3 789269 CY
4708531 2971173 247155.0 137212.1 10381130 CZ
4351269 3107740 318053.7 426690.0 82217837 DE
4321062 3653192 224577.6 177722.3 5475791 DK
5208731 4047869 171660.9 126292.2 1340935 EE
3180869 2021435 534870.4 440080.6 45283259 ES
template("", "maps/recmap1.png", basemap = FALSE, scale = FALSE)
plot.recmap(df, col = NA, border = col, lwd=4,  col.text = col)
dev.off()

cartog <- recmap(df)
template("", "maps/recmap2.png", basemap = FALSE, scale = FALSE)
plot(cartog[!cartog$name %in% c("IS","MT"),],  col = col, border = "white")
dev.off()

NB : IS et MT n’ont pu être placés. On les supprime à l’affichage.

cartogramR

library(cartogramR)

Le package cartogramR est déveoppé par Pierre-Andre Cornillon et Florent Demoraes de l’université de Rennes. Le package propose les méthodes flow based cartogram (Gastner and Newman 2004), fast flow based cartogram (Gastner, Seguy, and More 2018) et rubber band based cartogram (Dougenik, Chrisman, and Niemeyer 1985). La méthode flow based cartogram est basée sur go_cart.

Attention, cartogramR n’est pas/plus sur le CRAN. L’archive est néanmoins disponible et permet l’installation.

la fonction precartogramR() aide à choisir la taille de la grille de déformation (calcul un peu long). on choisit donc a minima la grille telle que le minimum d’intersections est superieur ou egal a un (ici un pas de grille de 256 peut faire l’affaire)

precarto <- precartogramR(communes, method = "GastnerSeguyMore")
summary(precarto)
##              L256      L512     L1024     L2048
## Min.      2.00000    5.0000   20.0000    82.000
## 1st Qu.  13.00000   52.0000  208.7500   837.000
## Median   19.00000   77.0000  306.0000  1226.500
## Mean     26.16992  104.6953  418.6816  1674.701
## 3rd Qu.  29.00000  115.2500  460.7500  1850.250
## Max.    405.00000 1613.0000 6469.0000 25884.000

la fonction cartogramR() parmet de calculer le cartogramme en tant que tel selon les 3 méthodes proposées : GastnerSeguyMore" (ou “gsm),”GastnerNewman" (ou “gn),”DougenikChrismanNiemeyer" (ou “dcn”). L’option L=256 permet de choisir la taille de la grille. L’option grid=TRUE (pour les méthodes “gsm” et “gn”) permet d’afficher la grille de déformation. L’option maxit=50 permet de définir le nombre d’itérations max (défaut = 50).

GastnerSeguy <- cartogramR(communes, count="pop2018", method="GastnerSeguyMore", options=list(L=256, grid=TRUE, maxit = 5))
## WARNING criterion: 26.455251 > Objective: 0.010000
##  Increase maxit or decrease Objective

La fonction make_layer() permet de récupérer la grille de calcul (mais aussi la grille d’origine, les centroides, etc)

grid <- make_layer(GastnerSeguy, type = c("final_graticule"))

Ensuite, il est très facile d’afficher le cartogram avec plot (via sf) ou comme précédemment, avec le package mapsf.

template("Gastner, Seguy & More", "maps/gastnerseguy.png", basemap = FALSE, scale = FALSE)
mf_map(GastnerSeguy$cartogram, col = col, border = "white", lwd = 1, add = TRUE)
mf_map(grid, col = NA, border = "#6b4266", lwd = 0.05, add = TRUE)
dev.off()

La fonction residuals() permet de calculer les erreurs liées à la déformation (erreur relative : taille finale / taille theorique * 100)

table = cbind(GastnerSeguy$initial_data[,c("id","pop2018")] %>% st_drop_geometry(),
      orig_area = GastnerSeguy$orig_area,
      final_area = GastnerSeguy$final_area,
      errors = residuals(GastnerSeguy, type = "relative error")*100
      )
knitr::kable(table[c(0:10),], row.names = F, digits = 1)
id pop2018 orig_area final_area errors
38001 6336 27223532 39376361.8 -0.2
38002 1025 16353553 5636619.1 -11.7
38003 1127 8082831 7027719.4 0.2
38004 1279 10663444 7921991.8 -0.5
38005 954 56515131 5540617.6 -6.7
38006 4062 34980514 24517963.6 -3.0
38008 28 4789504 95852.1 -45.0
38009 1009 5121670 6220553.5 -1.0
38010 686 13149817 4281634.4 0.3
38011 1073 8596299 6774560.2 1.4

Export au format sf

st_write(as.sf(GastnerSeguy),"files/GastnerSeguy.gpkg", delete_layer=TRUE)
# GastnerSeguy_sf = st_as_sf(table, geometry = GastnerSeguy$cartogram)
# st_write(GastnerSeguy_sf,"files/GastnerSeguy.gpkg", delete_layer=TRUE)

Variations

Dots Cartograms

La méthode des dots cartograms est une représentation cartographique à l’intersection des cartogrammes de Dorling et des cartes par points. Les premières cartes utilisant cette méthodes ont été réalisées sur les données du Covid (voir et voir). Article à paraitre. Voir aussi sur Observable

dotcartogram = function(x,var,itermax,onedot,radius){
crs = sf::st_crs(x)
coords <- st_coordinates(st_centroid(sf::st_geometry(x),of_largest_polygon = TRUE))
x <- x[c("id",var)] %>% st_drop_geometry()
x <- data.frame(x, coords)
colnames(x) <- c("id","v","x","y")
x$v <- round(x$v/onedot,0)
x <- x[x$v > 0,]
dots <- x[x$v == 1,c("x","y","v")]
rest <-  x[x$v  > 1,c("x","y","v")]

nb <- nrow(rest)
  for (i in 1:nb){
    new <- rest[i,]
    new$v <- 1
    for (j in 1:rest$v[i]){ dots <- rbind(dots,new)}
  }
  
dots$x <- jitter(dots$x)
dots$y <- jitter(dots$y)
dots$v <- radius
  
simulation <- circleRepelLayout(x = dots, xysizecols = 1:3,
                                  wrap = FALSE, sizetype = "radius",
                                  maxiter = itermax, weights =1)$layout

circles <- st_buffer(sf::st_as_sf(simulation, coords =c('x', 'y'),
                                    crs = crs), dist = radius) 
return(circles)
}
onedot = 1000
dc = dotcartogram(x = communes, var = "pop2018", itermax = 120,
                  onedot = onedot, radius = 600)
template(paste0("Un point = ",onedot," personnes"), "maps/dotcartogram.png", scale = FALSE)
mf_map(isere, col = "#CCCCCC", border = "white", lwd = 0.5, add = TRUE)
mf_map(dc, col = col, border = "#6b4266", lwd = 0.5, add = TRUE)
dev.off()

Anti Cartogramme

Hexbin…

A vous de jouer

Explorez les cartogrammes à l’échelle mondiale.

world <- st_read("data/world_countries.geojson", quiet = TRUE ) %>%
  st_transform( "+proj=bertin1953")

world <- world[world$ISO3 != "ATA",]
knitr::kable(world[c(0:10),], row.names = F, digits = 1)
ISO2 ISO3 ISONUM NAMEen NAMEfr UNRegion GrowthRate PopDensity PopTotal JamesBond geometry
GQ GNQ 226 Equatorial Guinea Guinée-Équatoriale Africa 2.9 30.1 845060 0 MULTIPOLYGON (((-514804.2 -…
TV TUV 798 Tuvalu Tuvalu Oceania 0.3 330.5 9916 0 MULTIPOLYGON (((11846963 52…
MV MDV 462 Maldives Maldives Asia 1.7 1212.2 363657 0 MULTIPOLYGON (((5515581 -20…
AD AND 20 Andorra Andorre Europe -1.9 149.9 70473 0 MULTIPOLYGON (((-1013908 15…
AG ATG 28 Antigua and Barbuda Antigua-et-Barbuda America 1.0 208.7 91818 0 MULTIPOLYGON (((-6420562 59…
AT AUT 40 Austria Autriche Europe 0.3 103.7 8544586 4 MULTIPOLYGON (((27563.98 73…
BE BEL 56 Belgium Belgique Europe 0.6 373.2 11299192 0 MULTIPOLYGON (((-621349.5 1…
BH BHR 48 Bahrain Bahren Asia 1.4 1812.2 1377237 0 MULTIPOLYGON (((2804520 -10…
BS BHS 44 Bahamas Bahamas America 1.2 38.8 388019 3 MULTIPOLYGON (((-6563129 25…
BB BRB 52 Barbados Barbade America 0.3 661.0 284215 0 MULTIPOLYGON (((-6530926 70…

Votre template cartographique

col = "#c291bc"
credits = "Vous, 2021"
theme = mf_theme(x = "default", bg = "#f0f0f0", tab = FALSE, 
                   pos = "center", line = 2, inner = FALSE, 
                   fg = col, mar = c(0,0, 2, 0),cex = 1.9)

templateworld = function(title, file){

  mf_export(
    world,
    export = "png",
    width = 1000,
    filename = file,
    res = 96,
    theme = theme, 
    expandBB = c(-.02,0,-.02,0)
)

  mf_title(title)

  mf_credits(
    txt = credits,
    pos = "bottomleft",
    col = "#1a2640",
    cex = 0.8,
    font = 3,
    bg = "#ffffff30"
  )
}

Votre carte de base

templateworld("Le monde en projection Bertin (thx Fil)", "maps/world.png")
mf_map(world, col =col, border = "white", lwd = 0.5, add = TRUE)
dev.off()

A vous de jouer…

Blibliographie

Dorling, D. 1996. “Area Cartograms: Their Use and Creation, Concepts and Techniques in Modern Geography.”
Dougenik, James A, Nicholas R Chrisman, and Duane R Niemeyer. 1985. “An Algorithm to Construct Continuous Area Cartograms.” The Professional Geographer 37 (1): 75–81.
Gastner, Michael T, and Mark EJ Newman. 2004. “Diffusion-Based Method for Producing Density-Equalizing Maps.” Proceedings of the National Academy of Sciences 101 (20): 7499–7504.
Gastner, Michael T, Vivien Seguy, and Pratyush More. 2018. “Fast Flow-Based Algorithm for Creating Density-Equalizing Map Projections.” Proceedings of the National Academy of Sciences 115 (10): E2156–64.
Heilmann, Roland, Daniel A Keim, Christian Panse, and Mike Sips. 2004. “Recmap: Rectangular Map Approximations.” In IEEE Symposium on Information Visualization, 33–40. IEEE.
Olson, Judy M. 1976. “Noncontiguous Area Cartograms.” The Professional Geographer 28 (4): 371–80.
Panse, Christian. 2016. “Rectangular Statistical Cartograms in r: The Recmap Package.” arXiv Preprint arXiv:1606.00464.